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ABSTRACT 


This  report  investigates  the  feasibility  of  a  numerical  solution  for  a  incom¬ 
pressible  viscous  flow  through  a  multiple-blade-row  turbomachinery.  The 
numerical  method  used  is  a  three-dimensional  incompressible  Navier-Stokes 
solver  based  on  flux-difference  splitting,  and  an  implicit  high-resolution  scheme. 

In  order  to  validate  the  numerical  predictions,  the  measurements  from  a  pre¬ 
swirl  High  Reynolds  Number  Pump  (HIREP)  experiment  are  selected  for 
the  purpose  of  comparisons.  Comparisons  between  predictions  and  measure¬ 
ments  are  presented.  The  comparisons  include  the  surface  pressure  distri¬ 
bution,  the  wake  behind  the  trailing  edge,  and  flow  visualization  on  the  hub 
and  blade  surfaces. 

ADMINISTRATIVE  INFORMATION 

This  investigation  was  authorized  by  the  Naval  Sea  Systems  Command  under  the 
NSSN  Propulsor  Development  Program,  in  accordance  with  Program  Element  63561N, 
Task  Area  52177002,  and  Work  Request  Number  N0002494WR1 032788.  This  work  was 
performed  under  DTMB  Work  Unit  5080-053. 

INTRODUCTION 

Due  to  the  recent  advancements  on  the  subject  of  the  Computational  Fluid  Mechan¬ 
ics  (CFD)  and  the  availability  of  high  speed  computing  machines,  various  Navier-Stokes 
solvers  are  used  to  simulate  relatively  complicated  flow  phenomenon  with  encouraging 
successes.  A  proper  application  of  an  accurate  and  reliable  numerical  procedure  can 
often  guide  the  designer  of  a  flow  machinery  to  identify  and  to  avoid  potential  diffi¬ 
culties,  and  help  to  select  the  most  favorable  and  innovative  configurations  for  model 
testings.  The  attractive  benefits  of  a  successful  simulation  are  reducing  the  cost  and  the 
cycling  time  of  the  design.  The  purpose  of  this  report  is  to  explore  the  feasibility  of  a 
numerical  solution  for  the  flow  through  a  complicated  and  practical  multiple-blade-row 
pump  configuration  and  to  evaluate  the  accuracy  of  the  solutions  by  comparing  them 
with  the  available  data. 

In  order  to  validate  the  numerical  predictions,  the  measurements  from  a  pre-swirl 
High  Reynolds  Number  Pump  (HIREP)  experiment  are  selected  for  the  purpose  of 
comparisons.  The  HIREP  experiment  was  conducted  at  the  the  Garfield  Thomas  water 
tunnel  at  Pennsylvania  State  University  under  the  sponsorship  of  the  Advanced  Re¬ 
search  Project  Agency  (ARPA).  One  of  the  purposes  of  the  experiment  is  to  provide  a 
data  base  for  numerical  simulations.  The  techniques  and  the  uncertainties  of  the  mea¬ 
surement  were  well  documented,^  hence,  the  reported  data  are  suitable  for  validating 
the  method  of  computational  fluid  dynamics  (CFD).  Because  both  the  separation  and 
a  significant  amount  of  the  secondary  flow  are  driven  by  viscous  effects  in  the  boundary 
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layer  at  the  end  wall  regions,  a  Reynolds  Averaged  Navier-Stokes  (RANS)  type  solver 
is  needed.  At  present,  a  three-dimensional  Navier-Stokes  solver  for  incompressible  flow 
based  on  an  implicit  high-resolution  finite-differencing  scheme  suggested  by  Harten,^ 
Yee  et  al.,^  and  Hartwich  et  al.^  is  employed  for  this  purpose.  The  solution  is  lim¬ 
ited  to  the  steady-state  case.  The  comparisons  between  predictions  and  measurements 
include  the  surface  pressure  distribution,  the  wake  behind  the  trailing  edge,  and  flow 
visualization  on  the  hub  and  blade  surfaces. 

Recently,  similar  viscous  simulations  were  performed  by  Dreyer  and  Zierke,^  and  by 
Lee  et  ah.®  To  obtain  a  numerical  solution,  the  former  used  a  cell-center  finite- volume 
technique  with  an  explicit  multistage  time  integration  procedure,^  and  the  latter  used 
a  pressure-based  implicit  relaxation  method.® 

DESCRIPTION  OF  EXPERIMENT 

The  experiment  described  next  is  one  of  many  that  were  conducted  at  the  HIREP 
facility.  The  data  obtained  from  this  experiment  are  compared  with  the  results  of 
the  numerical  simulation.  Only  the  relevant  informations  about  the  geometry  of  the 
configuration  and  uncertainties  of  the  measurements  are  given  in  detail. 

HIREP  is  a  pre-swirl  unit  with  13  IGV  blades  and  7  after  skew  rotor  blades.  The 
radii  of  the  hub  and  casing  are  constant  at  26.67  cm.  (10.5  in.)  and  53.34  cm.  (21  in.), 
respectively.  This  gives  an  annular  area  through  the  pump  of  0.67  m^  (7.22  ft^).  The 
inlet  guide  vanes  have  a  constant  chord  length  of  17.53  cm.  (6.90  in.)  and  a  solidity 
ranging  from  1.36  at  the  hub  to  0.68  at  the  tip.  The  rotor  blades  have  a  chord  length 
of  28.50  cm.  (11.22  in.)  and  solidity  of  1.19  at  the  hub,  and  a  chord  length  of  26.64 
cm.  (10.49  in.)  and  a  solidity  of  0.56  at  the  tip.  All  the  blades  have  fillets  at  the 
junctures.  The  gap  clearance  between  the  blade  tip  and  the  casing  is  0.130  in.  The 
distance  between  the  IGV  and  rotor  stack-up  lines  is  40.28  cm.  (15.86  in.)  Figure  1 
shows  a  set  of  IGV  and  rotor  blades  (looking  from  upstream  to  downstream).  The 
surface  definition  of  the  stator  blade  is  given  in  53  sections,  each  section  contains  4l 
nodes.  The  surface  definition  of  the  rotor  blade  is  given  in  53  sections,  each  section 
contains  36  nodes.  Grid  lines  shown  in  Fig.  1  were  used  to  create  a  smooth  surface  with 
a  B-spine  method.  Based  on  this  smooth  surface,  a  system  of  new  surface  grid  lines 
suitable  for  viscous  computation  was  then  generated. 

While  this  particular  experiment  was  conducted,  the  rotor  turned  at  its  design  point 
of  260  revolutions  per  minute  (rpm),  the  corresponding  speed  at  the  blade  tip  is  14.51 
m/s  (47.6  ft/s).  The  reference  velocity  used  as  a  monitor  velocity  during  the  entire 
experiment  is  10.63  m/s  (35.0  ft/s).  The  Reynolds  number  based  on  the  reference 
velocity  and  the  IGV  chord  length  is  2.3  million.  The  time  averaged  static  pressures 
on  both  the  IGV  and  the  rotor  blade  surfaces  were  measured  by  pressure  taps  which 
were  distributed  in  five  spanwise  rows  corresponding  to  nominal  locations  of  10,  30, 
50,  70,  and  90  percent  of  the  span.  The  overall  uncertainty  in  measurements  ranged 
from  413.64  pascal  (0.06  psi)  over  regions  of  minimal  streamwise  pressure  gradients  to 
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Fig.1.  Surface  definition  of  IGV  and  rotor  blades. 


896.22  pascal  (0.13  psi)  over  regions  of  large  streamwise  gradients.  A  radial  survey  of 
the  time  averaged  static  pressure,  stagnation  pressure  and  velocity  components  with  a 
five-hole  probe  were  conducted  at  37  percent  of  the  chord  length  upstream  of  the  IGV 
leading  edge  (defined  as  the  IGV  inlet  plane)  and  at  49.7  percent  chord  downstream 
of  the  IGV  trailing  edge  (defined  as  the  IGV  exit  plane).  At  a  given  circumferential 
location,  measurements  were  taken  at  50  different  radii.  In  order  to  obtain  better 
resolution  in  the  circumferential  direction,  a  five-hole  probe  rake  was  devised.  This 
rake  was  used  to  take  measurements  at  the  IGV  exit  plane  at  12  radial  locations.  At 
each  radial  location,  the  360  deg  survey  was  divided  into  966  increments  of  0.373  deg. 
The  percentage  uncertainty  for  the  five-hole  rake  is  2.6,  9.8,  and  9.9  percent  for  the 
axial  velocity,  the  radial  velocity,  and  the  tangential  velocity,  respectively. 

The  Laser  Doppler  Velocimetry  method  (LDV)  was  used  to  measure  the  axial  and 
tangential  velocity  components  both  upstream  and  downstream  of  the  rotor  blade  row. 
The  measurement  planes  were  located  at  26.7  percent  chord  axially  upstream  of  the 
rotor  tip  trailing  edge,  and  at  4.8,  21.4,  and  32.2  percent  chord  axially  downstream 
of  the  tip  of  the  rotor  blade  trailing  edge.  The  precision  error  in  mean  velocity  LDV 
measurements  ranged  from  OT  to  2.4  percent  depending  on  the  turbulence  level.  The 
measurement  of  radial  velocity  is  not  available. 

NUMERICAL  APPROXIMATION 

In  this  report,  the  centerline  of  the  hub  is  defined  as  the  X-axis  of  the  coordinate 
system,  with  the  positive  direction  being  from  upstream  to  downstream.  The  Z  axis 
is  chosen  to  lie  in  the  vertical  plane  and  is  positive  pointing  upwards.  The  Y  axis  is 
chosen  such  that  it  forms  a  right-handed  system  with  the  X  and  Z  axes.  The  positive 
tangential  velocity  points  in  the  direction  opposite  to  the  direction  of  the  rotating  rotor 
(HIREP  had  a  right-headed  rotor,  the  rotor  turned  in  clockwise  direction  while  looking 
from  downstream  to  upstream).  The  positive  radial  velocity  points  away  from  the  hub. 
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MOTION  RELATIVE  TO  A  ROTATING  REFERENCE  FRAME 

Relative  to  a  rotating  reference  frame,  Newton’s  equations  of  motion,  on  which  the 
Navier-Stokes  equations  are  derived,  are  not  directly  applicable.  In  analyses,  it  is  more 
convenient  to  specify  the  boundary  conditions  in  terms  of  the  rotating  frame  and  to 
modify  the  equations  of  motion  to  accommodate  the  transformation  to  a  non-inertial 
frame.  To  an  observer  fixed  with  the  rotating  frame,  the  expression  for  acceleration  is  no 
longer  an  invariant.  The  additional  terms  resulting  from  the  transformation  of  reference 
frames  appear  as  body  forces  and  can  be  accounted  for  by  introducing  centrifugal  and 
Coriolis  forces.  This  statement  can  be  expressed  mathematically  as  follows: 

where  u  is  the  velocity  vector,  12  is  the  angular  velocity  vector,  and  r  is  the  position 
vector.  The  subscripts  I  and  R  refer  to  inertial  and  rotating  frames  of  reference.  The 

term  is  the  actual  acceleration  that  a  fluid  particle  is  experiencing,  and 

is  the  acceleration  relative  to  the  rotating  frame;  12  x  (12  x  r),  and  212  x  ur  are  the 
centrifugal  and  Coriolis  forces,  respectively.  At  present,  the  flow  variables  in  the  rotor 
blade  row  are  computed  at  a  rotational  frame  of  reference. 

The  three-dimensional  incompressible  RANS  equations  based  on  primitive  variables 
are  formulated  in  a  boundary-fitted  curvilinear  coordinate  system.  Using  Chorin’s  arti¬ 
ficial  compressibility  formulation,®  the  incompressible  Navier-Stokes  equation  is  written 
in  conservation  law  form  for  three-dimensional  flow  with  the  reference  frame  rotating 
about  the  J^-axis  at  a  constant  angular  velocity  as, 

Qt  +  {E*  -  e:),  +  {F*  -  F:)y  +  (G*  -  G:),  -f  ^  =  O  .  (1) 

In  Eq.  1  the  dependent  variable  vector  Q  is  defined  as  Q  =  {p,  iz,  u,  w)^  and  the  inviscid 
flux  vectors  E*,F*,  and  G*;  the  viscous  shear  flux  vectors  E*,F*,  and  G*  ;  and  the 
body  forces,  H,  are  given  by 

E*  =  {I3u,u^  +  p,uv,uw)^ 

F*  =  {^v,uv,v^  +  p,vw)'^ 

G*  =  {/3w,uw,vw,w^  +  p)^ 

Ey  -  RG  (0,  Tjjj;,  Tjjy, 

E^  —  {^iTyxfTyy^'^yz) 

Gy  —  Rg  (0,  Tjxi ) ‘Tzz) 

H  =  (0, -12x^1/ +  Ditz;, -11x^2  -  0)^  . 
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The  coordinates  x,y,z  are  scaled  with  an  appropriate  characteristic  length  scale  L. 
The  Cartesian  velocity  components  u,v,w  are  nondimensionalized  with  respect  to  the 
free  stream  velocity,  Voo-  The  X  component  of  angular  velocity,  is  normalized 
with  Voo  and  T,  while  the  normalized  pressure  is  defined  as  p  =  (p  —  Pco)lpV^-  The 
kinematic  viscosity,  is  assumed  to  be  constant,  and  the  Reynolds  number  is  defined  as 
Re  =  V^Ljv.  The  artificial  compressibility  parameter,  j3,  monitors  the  error  associated 
with  the  addition  of  the  unsteady  pressure  term  dp/dt  in  the  continuity  equation  which 
is  needed  for  coupling  the  mass  and  momentum  equations  in  order  to  make  the  system 
hyperbolic. 

Equation  1  can  be  transferred  to  a  curvilinear,  body-fitted  coordinate  system  {C  ■> 
7}  )  through  a  coordinate  transformation  of  the  form 

C  =  C(a^,y,2),  (  =  i{x,y,z),  and  7/ =  p(a:,  y,  2)  . 

Thus,  Eq.  1  becomes, 

(Q/J),  +  (£-£„)<  +  (F  -  a)s  +  (G  -  G,),  +  (HjJ)  =  0  (2) 

with 

(E,F,Gf  =  {[T]  (E-jJ  ,F-IJ  ,G’jjf} 

and 

(E„F„G,f  =  [T\  (E-JJ  ,F;IJ  ,G-Jjf  , 

where 


[r|  = 


Cx  Cy  Cz 

^x  ^y 


L  Vx  Vy  Vz  J 


and  the  Jacobian  of  the  coordinate  transformation  is  given  by 


J  ^  =  det 


H  Vi 
H  Vi  H 

Vv  . 


The  Cartesian  derivatives  of  the  shear  fluxes  are  obtained  by  expanding  them  using 
chain  rule  expansions  in  the  C,  and  rj  directions. 

The  basic  operations  of  converting  the  set  of  differential  equations  to  a  system  of 
difference  equations  may  be  divided  into  temporal  differencing  and  spatial  differencing. 
The  procedure  is  described  in  the  following  sections. 
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TEMPORAL  DIFFERENCING 


Because  only  the  steady-state  solutions  are  of  interest,  a  first-order  accurate  Euler- 
implicit  time  differencing  scheme  is  used.  The  application  of  the  scheme  avoids  an 
overly  restrictive  time-step  size  when  highly  refined  grids  are  used  to  resolve  viscous 
effects.  In  addition,  a  spatially  variable  time  step  is  used  to  accelerate  convergence. 

Defining  computational  cells  with  their  centroids  at  I  =  6/A6  {  9  is  (,  (,ot  t]) 
and  their  cell  interfaces  at  /  ±  1/2,  the  backward  Euler  time  differencing  of  the  three- 
dimensional  conservation  form  is 

^  = -[Aj(£“+' - 

+  a,(G”«-g;+’)1-h/j  (3) 

where  At  is  the  time  step,  AQ"  =  and  Ai{E)  =  [(E);+i/2  —  {E)i_if2]/A(. 

Superscripts  denote  the  time  level  at  which  the  variables  are  evaluated. 


Linearizing  Eq.  .3  about  time  level  n,  we  obtain 


where 

RHS  =  -[A^(E"  -  F:)  +  A^(F"  -  F^)  +  A,(G’"  -  G^)]  -  ^/J 
and  where  I  is  the  identity  matrix. 

SPATIAL  DIFFERENCING 

The  three-dimensional  differential  operator  in  Eq.  4  is  first  split  into  three  indepen¬ 
dent  one-dimensional  operators.  The  spatial  differencing  of  the  inviscid  flux  in  each 
of  these  one- dimensional  operators  is  then  constructed  by  an  upwind  flux  differencing 
scheme  based  on  Roe’s  approximate  Riemann  solver  approach.^®  In  each  computational 
cell  the  differential  operator  is  linearized  around  an  average  state  such  that  the  flux 
difference  between  two  adjacent  cells  satisfies  certain  conservative  properties.  As  a  re¬ 
sult,  the  flux  at  an  interface  can  be  expressed  in  terms  of  the  directions  of  the  travelling 
waves.  Harten’s  high-resolution  total  variation  diminishing  (TVD)  technique^  is  then 
applied  to  enhance  the  accuracy  of  the  solution  to  a  higher  order  in  the  region  where 
its  variation  is  relatively  smooth.  The  undesirable  spurious  numerical  oscillations  as¬ 
sociated  with  higher  order  approximations  are  suppressed  by  a  limiter  suggested  by 
Wang  and  Richards. This  limiter  is  a  modified  Yee’s  version,^  and  is  more  suitable  for 
viscous  solutions  at  high  Reynolds  number.  This  approach  is  illustrated  in  the  following 
section. 
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Letting  the  flux  Jacobians  A,  B,  and  C  be  deflned  as  follows: 

^  "  IQ'  ^  ^  dQ'^  ""  dQ' 

discretize  the  inviscid  and  viscous  fluxes  according  to  the  upwind  differencing  scheme 
and  to  the  central  differencing  scheme,  respectively,  in  C>  ti  ^nd  tj  coordinate  directions 
independently  and  then  assemble  them  together.  Equation  4  becomes 

lil/JM)  -  (/I- +  A,+i  +  (A*  + 

-  (B-  +  K);^jA,-+i+(B+  +  y),.,A,_x 

-  (C-  +  Z)x+i Ax+j  +  (C*  +  Z)x_iAx_xl"AO“ 

=  -RESiQ")  (5) 

where  i,j,  and  k  are  spatial  indices  associated  with  the  r/,  and  (  coordinate 
directions,  and  ^=*=,5^,  and  are  flux  matrices  split  from  the  flux  Jacobians  A,B, 
and  C  according  to  the  signs  of  their  eigenvalues.  The  residual  RES{Q'^)  is  evaluated 
with  a  TVD  technique;  the  discretization  is  accurate  to  the  third  order.  Conventional 
second-order  central  differencing  is  applied  to  obtain  the  viscous  flux  matrices  X,  Y, 
and  Z.  Equation  5  is  solved  by  an  implicit  hybrid  algorithm,  where  a  symmetric 
planar  Gauss- Seidel  relaxation  is  applied  in  the  streamwise  direction  C  in  combination 
with  approximate  factorization  in  the  remaining  two  coordinate  directions  (  and  t)  to 
avoid  the  At^  spatial  splitting  error  incurred  in  fully  three-dimensional  approximate 
factorization  methods.  This  scheme  is  unconditionally  stable  for  linear  systems  and 
can  be  vectorized  completely.  As  a  result  of  upwind-differencing,  the  coefficient  matrix 
of  the  system  becomes  diagonally  dominant.  In  addition,  the  necessity  of  adding  and 
tunning  a  numerical  dissipation  term  for  stability  reasons,  as  in  some  schemes  with 
central  differencing,  is  alleviated. 

TURBULENCE  MODEL 

The  Baldwin- Lomax  algebraic  turbulence  model  (BLM)  is  modified  to  simulate  the 
turbulent  effect  of  the  flow.  The  BLM  is  applied  to  the  blade- to-blade  (^)  and  spanwise 
(rj)  directions  independently.  The  mixing  length  I  for  the  inner  layer  is  calculated  based 
on  the  distances  from  the  end- wall  {s^)  and  from  the  blade  surface  (s^): 

"t~  Sjj  “H  “b 

As  suggested  by  Chima  et  al.^^  and  Dorney  and  Davis, the  outer  and  the  inner  tur¬ 
bulence  viscosities  are  combined  by  using  Granville’s  blending  function, that  is,  the 
combined  viscoscity  //  is  set  to 

Uotanh— 

/^o 
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where  /Zq  and  //,•  are  the  outer  and  inner  viscosities,  respectively.  The  turbulence  veloc¬ 
ity  in  the  corner  flow  region  is  computed  based  on  a  blending  function  by  Vatsa  and 
Wedan^®: 


_  Mg  T  ■Sg  firi 

2  I  2 

In  addition,  the  location  of  the  wake  center  line  is  searched  before  the  application  of 
the  BLM. 


GRID  TOPOLOGIES  AND  BOUNDARY  CONDITIONS 

Grid  topologies  used  to  compute  the  IGV  and  rotor  passage  flow  are  similar.  Essen¬ 
tially,  the  system  includes  three  main  blocks.  There  is  a  block  of  0-type  grid  around 
the  blade;  it  is  then  connected  upstream  and  downstream  to  two  blocks  of  H-type  grid 
that  are  extended  to  upstream  boundary  and  downstream  boundary  planes,  respec¬ 
tively.  Through  the  overlapping  cells  at  the  connecting  region,  the  solutions  between 
the  blocks  communicate.  In  case  of  rotor  passage,  an  additional  block  of  H-type  grid 
is  used  to  cover  the  tip  gap.  For  the  0-type  grid,  there  are  227  nodes  around  the 
blade  surface,  33  nodes  in  the  direction  normal  to  the  blade  surface  and  95  nodes  from 
hub  to  casing.  The  upstream  H-type  grid  contains  21  nodes  in  the  axial  direction,  39 
nodes  in  the  circumferential  direction,  and  95  nodes  in  the  radial  direction.  The  size 
of  the  downstream  H-type  grid  is  21  x  77  x  95  in  the  axial,  circumferential,  and  radial 
directions,  respectively.  The  total  number  of  the  nodes  for  each  blade  passage  is  about 
943,000. 

The  flllets  at  the  blade  hub/casing  junctures  are  included  in  the  gridding  process. 
The  grid  spacings  near  the  solid  surface  are  selected  such  that  the  the  first  y'*'  coordi¬ 
nate  off  the  surface  is  about  1  or  less.  Figure  2  shows  the  surface  grids  in  blade  passages 
with  every  other  grid  in  each  coordinate  direction  removed.  The  periodicity  conditions 
are  applied  along  the  bounding  surfaces  in  the  circumferential  direction.  At  the  up¬ 
stream  boundary  plane,  the  velocity  components  are  specified  and  the  static  pressure  is 
extrapolated  from  the  interior.  At  the  downstream  boundary  plane,  the  static  pressure 
is  defined  and  the  velocity  components  are  extrapolated  from  the  interior.  On  the  solid 
surfaces,  all  components  of  velocity  and  the  normal  gradient  of  static  pressure  vanish.  In 
addition,  in  the  rotor  passage,  the  tangential  velocity  on  the  surface  of  the  nonrotating 
portion  of  the  hub  is  prescribed  according  to  the  rotating  speed  of  the  rotor. 

SOLUTION  PROCEDURE 

In  order  to  conserve  the  computational  effort,  the  plenum  region  of  the  water  tunnel 
is  not  modeled,  only  the  region  where  the  hub  and  casing  are  cylindrical  is  covered. 
The  domain  extends  approximately  from  1.5  IGV  chord  lengths  upstream  of  the  IGV 
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Fig.  2a.  IGV  passage. 


Fig.  2.  Rotor  passage. 


Fig.  2.  Surface  grids  in  blade  passages. 


leading  edge  to  approximately  2.5  IGV  chord  lengths  downstream  of  the  rotor  blade 
trailing  edge.  All  the  measurements  fall  within  this  domain.  Because  only  steady  state 
solutions  are  of  interest,  a  simplified  iterative  procedure  is  proposed.  The  precedure  is 
described  as  follows: 

1.  Establish  a  computational  domain  including  the  IGV  blade- to-blade  passage. 
The  grid  system  extended  from  1.5  IGV  chord  lengths  upstream  of  the  blade  leading 
edge  to  about  2.5  IGV  chord  lengths  downstream  of  the  IGV  trailing  edge.  The  inlet 
velocity  distribution  is  assumed  to  be  axisymmetrical  and  is  similar  to  the  l/7‘*  power 
boundary  layer  profile.  It  is  then  scaled  such  that  the  mass  flow  rate  agrees  with  the 
measurement  taken  at  37  percent  IGV  chord  upstream  of  the  blade  leading  edge.  The 
rotor  effect  can  be  approximated  by  distributing  the  body  forces  to  the  computational 
cells  located  within  the  rotor  swept  volume.  The  magnitude  of  the  body  forces  is 
determined  analytically  based  on  the  rotor  characteristics.  Use  the  grid  system  and  the 
boundary  conditions  described  earlier  to  obtain  a  converged  solution. 

A  plane  at  49.7  percent  chord  behind  the  trailing  edge  is  defined  as  the  stator  exit 
plane.  Circumferentially  averaged  solutions  at  that  plane  are  then  used  as  input  for  the 
rotor  inlet  plane. 

2.  Establish  a  computational  domain  including  the  rotor  blade-to-blade  passage. 
The  grid  system  extends  from  0.75  IGV  chord  lengths  upstream  of  the  blade  leading 
edge  to  2.5  IGV  chord  lengths  downstream  of  the  trailing  edge.  The  axisymmetrical 
velocity  profiles  at  the  inlet  plane  are  defined  from  the  circumferentially  averaged  solu¬ 
tion  obtained  at  the  exit  plane  in  step  1.  The  solutions  in  the  passage  are  then  obtained 
by  solving  the  equations  of  motion  in  the  rotating  reference  frame. 

Use  the  newly  obtained  circumferentially  averaged  static  pressure  at  the  inlet  plane 
from  step  2  as  the  pressure  boundary  condition  at  the  exit  plane  for  step  1  and  repeat 
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the  process.  At  present,  the  effort  of  performing  a  second  iteration  is  minimum,  because 
the  interaction  between  the  blade  rows  is  small. 

Computations  were  performed  on  a  CONVEX  3080  mini  supercomputer.  For  930,000 
nodes,  each  iteration  required  2  minutes  single  processor  CPU  time.  For  each  blade  pas¬ 
sage,  it  required  500  iterations  to  rearch  converaged  solutions  (L2  norm  of  the  residual 
drops  by  about  three  order  of  magnitude).  To  obtain  above  results,  700  iterations  were 
carried  out. 


RESULTS 

The  comparisons  between  the  numerical  results  and  the  measurements  are  presented 
in  the  following  sections. 

IGV  PASSAGE 

Due  to  the  interaction  of  boundary  layers  from  blades  and  wall  surfaces  (such  as 
the  hub  and  casing),  the  secondary  flow  is  generated  as  the  flow  passes  through  the 
blade  row  passage.  Secondary  flow  often  appears  in  the  form  of  vortices;  common 
ones  are  horseshoe  vortex,  trailing  edge  vortex,  passage  vortex,  corner  vortex,  and  tip 
gap  vortex.  In  this  report,  the  secondary  flow  at  a  given  plane  is  defined  in  the  same 
manner  as  described  by  Zierke  et  al.^  The  local  area-weighted  circumferential  mean  at 
each  spanwise  location  was  computed  and  then  was  substracted  from  the  cross  flow 
components  at  the  same  spanwise  location.  Figure  3a  shows  that  some  of  the  vortices 
mentioned  previously  can  be  identified  clearly  from  the  measurements  taken  at  the 
HIREP  experiment  at  49.7  percent  chord  downstream  of  the  trailing  edge  of  the  IGV. 
The  computational  result  at  the  same  location  is  shown  in  Fig.  3b.  As  the  flow  passes  a 
blade  section,  the  streamlines  on  the  suction  side  of  the  blade  tend  to  flow  inward  toward 
the  location  where  the  maxinum  lift  occurs,  while  the  streamlines  on  the  pressure  side 
of  the  blade  tend  to  flow  outward  from  the  same  location.  As  a  result,  a  surface  of 
discontinuity  represented  by  a  vortex  sheet  is  formed  as  the  flow  leaves  the  trailing  edge 
of  the  blade. 

The  IGV  blade  is  designed  such  that  the  maximum  lift  occurs  at  42  percent  span. 
From  the  formation  of  the  trailing  edge  vortex.  Fig.  3b  indicates  that  the  computa¬ 
tion  predicts  this  location  quite  well.  However,  physically,  the  trailing  vortex  sheet  is 
unstable.  At  a  pertubation  the  sheet  will  roll  up  into  a  pair  of  vortex  tubes  rotating 
in  opposite  directions.  Measurements  shown  in  Fig.  3a  suggest  that  such  a  process  is 
taking  place. 

In  the  core  region  of  the  blade  passage,  a  pressure  gradient  towards  the  suction 
surface  is  established  to  balance  the  circumferential  acceleration.  As  a  first  order  ap¬ 
proximation,  it  can  be  assumed  that  the  same  pressure  gradient  exists  throughout  the 
boundary  layers  near  the  hub  and  casing  surfaces.  For  a  given  pressure  gradient,  the 
smaller  the  velocity,  the  smaller  the  radius  that  the  particle  will  follow;  this  is  the  case 
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Fig.  3a.  Measurements  with  five-hole  probe. 


Fig.  3b.  Computation. 


Fig.  3.  Secondary  velocity  vectors  49.7  percent  chord  axially  downstream 

of  the  IGV  trailing  edge. 


in  the  wall  boundary  layers.  Therefore,  near  the  walls,  the  slower  moving  flow  is  forced 
toward  the  suction  side  from  the  pressure  side.  This  mechanism  generates  two  passage 
vortices  of  opposite  sign  between  the  blades  that  can  be  found  in  Fig.  3a  and  Fig.  3b. 
As  discussed  in  Zierke  et  al.,^  the  trailing  vortex  sheets  and  the  passage  vortices  coexist. 
One  is  not  actually  the  cause  of  the  other.  For  the  present  case,  the  strengths  of  the 
passage  vortices  are  weak,  as  can  be  indicated  by  the  flow  visualizations  on  the  hub 
and  casing  surfaces.  Boundary  layers  at  hub/casing  surfaces  separate  in  front  of  the 
leading  edge  of  the  blades,  and  horseshoe  vortices  are  generated.  The  pressure  side  lag 
of  the  horseshoe  vortex  has  the  same  sign  as  the  passage  vortex  and  will  merge  into  the 
passage  vortex.  Figures  2  and  3  present  no  evidence  of  the  horseshoe  vortices  at  49.7 
percent  chord  downstream  from  the  trailing  edge.  Perhaps,  their  strengths  are  greatly 
reduced  through  the  interaction  with  other  vortices. 

Figure  4  shows  the  comparison  of  the  measured  and  computed  IGV  wakes  at  sev¬ 
eral  spanwise  locations.  The  measuring  plane  is  49.7  percent  chord  axially  downstream 
of  the  IGV  trailing  edge.  The  streamwise  velocity  V*  is  defined  as  the  vector  sum  of 
the  velocities  in  the  axial  (14)  and  tangential  directions  (Vi).  The  experimental  data 
presented  are  passage-averaged.  Error  bars  are  computed  and  reported  in  Zierke  et  al.^ 
Near  the  hub  at  4.8  percent  span  and  9.5  percent  span  the  wakes  are  wider,  deeper, 
and  asymmetric  due  to  the  boundary  layer  interaction  at  the  blade/hub  region.  Com¬ 
putations  display  the  same  characteristics,  and  the  agreement  with  the  measurements 
is  very  good.  Farther  away  from  the  hub  surface,  the  wakes  become  shallower  and  more 
symmetric.  Computations  show  the  same  trend  but  with  somewhat  deeper  deficits.  Fig¬ 
ure  5  shows  the  circumferentially- averaged  velocity  components  at  49.7  percent  chord 
axially  downstream  of  the  IGV  trailing  edge.  Two  sets  of  measurement  are  shown  in  the 
figure,  one  was  obtained  with  radial  probe  and  the  other  with  wake  rake.  Uncertainty 
analysis  by  Zierke  et  al.^  shows  that  the  wake  rake  data  are  more  reliable. 
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Fig  4.  Passage-averaged  IGV  wakes  49.7  percent  chord  axially  downstream 

of  the  IGV  trailing  edge. 


Hollow  Symbols  -  Measurements  from  Radial  Probe 
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Lines  -  Computations 


Fig.  5.  Circumferentially  averaged  velocities  49.7  percent  chord  axially  downstream 

of  the  IGV  trailing  edge. 
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Figure  6  shows  the  static  pressure  distribution  on  the  IGV  blade  surface.  The  pres¬ 
sure  coefficient  is  normalized  by  a  dynamic  pressure  based  on  blade  tip  speed  as  defined 
by  Zierke  et  al.^  At  10  percent  span,  near  the  leading  edge,  the  measured  pressures 
at  both  suction  and  pressure  side  are  lower  than  those  computed.  The  same  discrep¬ 
ancy  was  found  in  the  lifting  surface  calculation^  and  in  Dreyer  and  Zierke’s  viscous  flow 
computation.^  Computations  agree  well  with  measurements  at  other  spanwise  locations. 

Figure  7  shows  the  flow  visualization  on  the  hub  surface.  The  locations  of  the  saddle 
points  at  the  hub  surface  near  the  leading  and  trailing  edges  are  predicted  correctly. 
Figure  8a  shows  the  flow  visualization  on  the  blade  pressure  side  surface.  Both  com¬ 
putation  and  experiment  show  no  flow  separation  on  the  surface.  Figure  8b  shows  the 
flow  visualization  on  the  blade  suction  side  surface.  Experiments  indicate  that  flow  sep¬ 
arates  upstream  of  the  trailing  edge  with  corner  vortices.  The  predicted  separation  line 
is  slightly  downstream  of  the  obsverved  line  except  near  the  casing  where  the  predicted 
separation  is  further  upstream  than  that  observed. 

ROTOR  PASSAGE 

Figure  9  shows  the  comparison  of  the  measured  and  computed  rotor  wakes  in  differ¬ 
ent  spanwise  locations.  The  measuring  plane  is  32.3  percent  chord  axially  downstream 
of  the  rotor  tip  trailing  edge.  At  each  spanwise  location,  the  axial  (\4)  and  tangential 
(Vi)  velocity  components  are  measured  with  the  LDV  method  around  360  deg  at  0.5  deg 
intervals.  The  spatial  repeatability  of  the  measurements  among  all  seven  blade  passages 
is  excellent;  only  two  cycles  of  representative  data  are  shown  in  Fig.  9.  The  agreement 
between  the  measurement  and  the  computation  is  very  good  from  20  percent  span  to 
near  80  percent  span.  Note  that  the  maximum  lift  on  the  blade  is  designed  to  occur  at 
about  22  percent  of  span.  At  locations  closer  to  the  hub  surface,  computed  wakes  have 
larger  than  measured  deficits.  At  locations  near  the  casing  end-wall,  the  velocity  peak 
at  mid-pitch  is  not  predicted  very  well. 

Figure  10  shows  the  static  pressure  distribution  on  the  rotor  surface.  Figure  11 
shows  the  computed  secondary  velocity  vectors  (Fig.  11a)  and  streamlines  (Fig.  llb)at 
a  plane  37.8  percent  chord  axially  downstream  of  the  rotor  tip  trailing  edge  (looking  from 
downstream  to  upstream).  The  blade  wake  can  be  identified  by  the  locations  where  the 
vectors  have  large  negative  (in  the  direction  of  blade  rotation)  tangential  components. 
Due  to  the  relative  motion  between  the  blade  and  the  end-wall,  the  blade  scrapes  up  the 
end-wall  boundary  layer  in  front  of  it  as  it  advances.  As  a  result,  a  scraping  vortex  is 
formed.^®  In  Fig.  11,  the  trace  of  the  scraping  vortex  can  be  located  ahead  of  the  blade 
wake.  Trailing  behind  the  blade  wake  near  the  end-wall  is  the  tip  leakage  vortex.  The 
passage  vortices  between  the  blades  and  the  horseshoe  vortex  system  at  the  blade/hub 
juncture  can  also  be  identified.  Figure  12  shows  the  flow  visualization  on  hub  surface. 
The  location  of  the  saddle  point  ahead  of  the  leading  edge  is  predicted  correctly.  The 
predicted  location  of  the  saddle  point  near  the  trailing  edge  is  further  downstream 
than  observed.  Figure  13  shows  the  flow  visualization  on  the  blade  suction  surface. 
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Fig.  9.  Circumferential  variation  of  velocity  components  32.2  percent  chord  axially  downstream 

of  the  rotor  tip  trailing  edge. 
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Fig.  10.  Static  pressure  distribution  on  the  rotor  surface. 


Fig.  11.  Predicted  secondary  veiocity  at  37.8  percent  chord  axialiy  downstream 

of  the  rotor  tip  traiiing  edge. 
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It  indicates  that  trailing  edge  separation  vortex  moves  radially  up  the  trailing  edge. 
The  predicted  separation  line  agrees  well  with  the  observed  line.  Figure  14  shows  the 
particle  traces  of  flow  through  the  tip  gap  from  the  pressure  side  to  the  suction  side 
of  the  blade.  Vorticity  is  formed  and  shedded  along  the  suction  surface.  A  similar 
phenomenon  was  observed  in  an  experiment  by  FarelF^  and  predicted  by  Dreyer  and 
Zierke.^ 


Fig.  14.  Particle  traces  of  the  flow  through  rotor  tip  gap. 

CONCLUSIONS  AND  RECOMMENDATIONS 

An  implicit  high-resolution  finite-differencing  scheme  to  obtain  a  steady-state  so¬ 
lution  for  flow  passing  a  multiple-blade-row  turbomachinery  is  shown.  In  the  scheme 
flows  through  stator  and  rotor  passages  were  solved  separately,  and  then  were  coupled 
through  iteration  processes. 

The  surface  flow  visualizations  indicate  that  the  flow  separations  on  the  hub  and 
blade  surfaces  are  predicted  correctly,  although  the  details  of  the  surface  imprint  of 
the  corner  vortices  near  the  junctures  are  not  captured.  The  visuahzations  further 
suggest  that  the  flow  near  the  end-wall  surfaces  turns  from  the  pressure  side  of  the 
blades  toward  the  suction  side  of  the  blades.  However,  the  turning  motion,  that  can 
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be  associated  with  the  formation  of  a  passage  vortex,  is  relatively  mild  in  comparison 
with  those  occurring  in  the  turbine  passages.  Because  of  the  relative  weakness  of  the 
passage  vortex,  the  predicted  secondary  flow  pattern  behind  the  IGV  agrees  well  with 
the  inviscid  solution  presented  by  Zierke  et  al.^  In  the  rotor  passage,  because  of  the 
relative  motion  between  the  rotor  blade  and  the  casing,  the  tip  gap  vortex  behind  the 
blade  and  the  scraping  vortex  in  front  of  the  blade  are  predicted. 

CONCLUSIONS 

Comparisons  between  the  computed  and  the  measured  static  pressures  on  the  IGV 
and  the  rotor  blade  surfaces  show  that  the  loading  distributions  are  predicted  reason¬ 
ably  well.  The  wakes  behind  the  IGV  are  predicted  reasonably  well  even  near  the  hub 
juncture.  Overall,  the  predicted  wake  deficits  are  somewhat  deeper  than  those  mea¬ 
sured.  The  predicted  and  measured  rotor  wakes  disagree  near  the  hub  and  casing,  but 
better  agreements  exists  at  mid-span  region. 

RECOMMENDATIONS 

The  computations  for  rotor  passage  were  carried  out  in  a  rotating  reference  frame. 
In  order  to  improve  the  prediction,  the  effects  of  rotation  to  the  turbulence  modeling 
require  some  special  considerations.  Some  useful  suggestions  can  be  found  in  a  report 
by  Gorski.^® 
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